home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
The World of Computer Software
/
The World of Computer Software.iso
/
yamp16.zip
/
VIRTOP.CPP
< prev
next >
Wrap
C/C++ Source or Header
|
1993-01-11
|
41KB
|
1,647 lines
/*
YAMP - Yet Another Matrix Program
Version: 1.6
Author: Mark Von Tress, Ph.D.
Date: 01/11/93
Copyright(c) Mark Von Tress 1993
Portions of this code are (c) 1991 by Allen I. Holub and are used by
permission
DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
FROM USE OF THIS PROGRAM.
*/
#include "virt.h"
/////////////////// matrix functions ---- non-member functions
int Setwid(int w)
{
static int wid = 6;
wid = (w < 0) ? wid : w;
return wid;
}
int Setdec(int d)
{
static int dec = 3;
dec = (d < 0) ? dec : d;
if (d >= Getwid()) dec = Getwid() - 1;
dec = (dec <= 0) ? 0 : dec;
return dec;
}
int Getwid(void)
{
static int wid = 6;
wid = Setwid(- 1);
return wid;
}
int Getdec(void)
{
static int dec = 3;
dec = Setdec(- 1);
return dec;
}
VMatrix &Reada(char *fid) /* this reads an ascii matrix */
{
unsigned int chh[MAXCHARS], *chc;
int i = 0, j = 0, x1 = 0, x2 = 0;
FILE *file;
char mname[MAXCHARS];
float x = 0;
file = fopen(fid, "r");
if (!file) Dispatch->Nrerror("READA: error while opening file");
for (i = 0; i < 80; i++) chh[i] = 0;
i = 0;
do { /* read until first carriage return */
j = fgetc(file);
chh[i] = j;
i++;
} while (j != ((int) 0x0A) && !ferror(file) && (i < 80));
if (ferror(file)) Dispatch->Nrerror(
"READA: error while reading matrix name" );
chh[i] = (int) '\0'; /* null terminate the string */
if (chh[0] != (int) 0x0A) {
strcpy(mname, ((const char *) chh));
chc = chh;
for (j = 1; j < i - 1; j++) strcat(mname, ((const char *) (++ chc)));
}
else {
strcpy(mname, " ");
}
fscanf(file, "%d%d", &x1, &x2);
VMatrix *mat = new VMatrix(mname, x1, x2);
for (i = 1; i <= x1; i++) {
for (j = 1; j <= x2; j++) {
fscanf(file, "%f", &x);
mat->M(i, j) = (double) x;
}
}
fclose(file);
Dispatch->Push(*mat);
delete mat;
return Dispatch->ReturnMat();
} /* reada */
VMatrix &Readb(char *fid) /* this reads an binary matrix */
{
int i = 0, j = 0, x1 = 0, x2 = 0;
FILE *file;
char name[MAXCHARS];
double x = 0;
file = fopen(fid, "rb");
if (!file) Dispatch->Nrerror("READB: error while opening file");
fread(&i, sizeof (int), 1, file);
fgets(name, (i + 1), file);
fread(&x1, sizeof (int), 1, file);
fread(&x2, sizeof (int), 1, file);
VMatrix *mat = new VMatrix(name, x1, x2);
for (i = 1; i <= x1; i++) {
for (j = 1; j <= x2; j++) {
fread(&x, sizeof (double), 1, file);
mat->M(i, j) = x;
}
}
fclose(file);
Dispatch->Push(*mat);
delete mat;
return Dispatch->ReturnMat();
} /* readb */
void Writea(char *fid, VMatrix &mat) /* write a matrix in an ascii file */
{
int i, j;
FILE *file;
file = fopen(fid, "w");
if (!file) Dispatch->Nrerror("WRITEA: unable to open file");
mat.Garbage("Writea");
fprintf(file, "%s\n", mat.StringAddress());
fprintf(file, "%d %d \n", mat.r, mat.c);
for (i = 1; i <= mat.r; i++) {
for (j = 1; j <= mat.c; j++)
fprintf(file, "%lf ", mat.m(i, j));
fprintf(file, "\n");
}
fclose(file);
} /* writea */
void heapsort(VMatrix &a)
{
int n = a.r;
int s0, s1, s2, s3, s4, s5, s4m1, s5m1;
double t1, ts;
s0 = s1 = n; /* Shell-Metzger sort, PC-World, May 1986 */
s1 /= 2;
while (s1 > 0) {
s2 = s0 - s1;
for (s3 = 1; s3 <= s2; s3++) {
s4 = s3;
while (s4 > 0) {
s5 = s4 + s1;
s4m1 = s4;
s5m1 = s5;
if (a.m(s4m1, 1) > a.m(s5m1, 1)) {
t1 = a.m(s4m1, 1);
a.M(s4m1, 1) = a.m(s5m1, 1);
a.M(s5m1, 1) = t1;
ts = a.m(s4m1, 2);
a.M(s4m1, 2) = a.m(s5m1, 2);
a.M(s5m1, 2) = ts;
s4 = s4 - s1;
}
else {
s4 = 0;
}
}
}
s1 /= 2;
}
}
VMatrix &MSort(VMatrix &a, int col, int order)
{
a.Garbage();
VMatrix *sorted = new VMatrix("sorted(", a.r, a.c);
if (!order) {
// sort column col, then reorder other rows
VMatrix *temp = new VMatrix("temp", a.r, 2);
for (int i = 1; i <= a.r; i++) {
temp->M(i, 1) = a.m(i, col);
temp->M(i, 2) = (double) i;
}
heapsort(*temp);
for (i = 1; i <= a.r; i++) {
for (int j = 1; j <= a.c; j++)
sorted->M(i, j) = a.m(((int) temp->m(i, 2)), j);
}
delete temp;
}
else {
// sort row col, then reorder other columns
VMatrix *temp = new VMatrix("temp", a.c, 2);
for (int i = 1; i <= a.c; i++) {
temp->M(i, 1) = a.m(col, i);
temp->M(i, 2) = (double) i;
}
heapsort(*temp);
for (i = 1; i <= a.c; i++) {
for (int j = 1; j <= a.r; j++)
sorted->M(j, i) = a.m(j, ((int) temp->m(i, 2)));
}
delete temp;
}
strtype name = sorted->Getname(*sorted);
name = name + a.Getname(a) + ")";
sorted->Nameit(name);
Dispatch->Push(*sorted);
delete sorted;
return Dispatch->ReturnMat();
}
VMatrix &Submat(VMatrix &a, int r, int c, int lr, int lc)
/* extract a submatrix of a into b*/
{
a.Garbage("Submat");
if ((r > a.r) || (c > a.c))
a.Nrerror("SUBMAT: index out of range");
VMatrix *b = new VMatrix("b", r - lr + 1, c - lc + 1);
strtype temp = a.Getname(a);
b->Nameit(temp);
int r2 = r - lr + 1;
int c2 = c - lc + 1;
for (int i = 1; i <= r2; i++) {
for (int j = 1; j <= c2; j++)
b-> M(i, j) = a.m(lr - 1+ i, lc - 1+ j);
}
Dispatch->Push(*b);
delete b;
return Dispatch->ReturnMat();
} /* submat */
VMatrix &Ch(VMatrix &a, VMatrix &b) /* horizontal concatination */
{
a.Garbage("Ch");
b.Garbage("Ch");
int adim = a.c;
int bdim = b.c;
int k = adim + bdim;
if (a.r != b.r)
a.Nrerror("CH: matrices have different number of rows");
VMatrix *c = new VMatrix("(", a.r, k);
if (!c) a.Nrerror("CH: failed to create temp storage");
strtype mname = c->Getname(*c) + a.Getname(a) + "||" + b.Getname(b) +")";
c->Nameit(mname);
for (int i = 1; i <= a.r; i++) {
for (int j = 1; j <= a.c; j++)
c->M(i, j) = a.m(i, j);
}
int ii;
for (i = 1, ii = 1; i <= b.r; i++, ii++) {
for (int j = 1, jj = 1; j <= b.c; j++, jj++)
c->M(ii, jj + a.c) = b.m(i,j);
}
Dispatch->Push(*c);
delete c;
return Dispatch->ReturnMat();
} /* ch */
VMatrix &Cv(VMatrix &a, VMatrix &b) /* vertical concatination */
{
a.Garbage("Cv");
b.Garbage("Cv");
int adim = a.r;
int bdim = b.r;
int k = adim + bdim;
if (a.c != b.c)
a.Nrerror("CV: matrices have different number of columns");
VMatrix *c = new VMatrix("(", k, a.c);
strtype mname = c->Getname(*c) + a.Getname(a) + "//" + b.Getname(b) +")";
c->Nameit(mname);
for (int i = 1; i <= a.r; i++) {
for (int j = 1; j <= a.c; j++)
c->M(i, j) = a.m(i, j);
}
int ii;
for (i = 1, ii = 1; i <= b.r; i++, ii++) {
for (int j = 1, jj = 1; j <= b.c; j++, jj++)
c->M(ii + a.r, jj) = b.m(i,j);
}
Dispatch->Push(*c);
delete c;
return Dispatch->ReturnMat();
} /* cv */
////////////////// math funtions
VMatrix& Tran(VMatrix &ROp)
{
ROp.Garbage("Tran");
VMatrix *temp = new VMatrix("t", ROp.c, ROp.r);
strtype newname = temp->Getname(ROp);
for (int i = 1; i <= ROp.r; i++) {
for (int j = 1; j <= ROp.c; j++) temp->M(j, i) = ROp.m(i, j);
}
newname = newname + "'";
temp->Nameit(newname);
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
VMatrix& Mexp(VMatrix &ROp)
{
// take exponent of all elements
ROp.Garbage("Mexp");
VMatrix *temp = new VMatrix("exp(", ROp.r, ROp.c);
for (int i = 1; i <= ROp.r; i++) {
for (int j = 1; j <= ROp.c; j++) temp->M(i, j) = exp(ROp.m(i, j));
}
strtype newname = temp->Getname(*temp) + ROp.Getname(ROp) + ")";
temp->Nameit(newname);
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
VMatrix& Mabs(VMatrix &ROp)
{
// take absolute values of all elements
ROp.Garbage("Mabs");
VMatrix *temp = new VMatrix("abs(", ROp.r, ROp.c);
for (int i = 1; i <= ROp.r; i++) {
for (int j = 1; j <= ROp.c; j++) temp->M(i, j) = fabs(ROp.m(i, j));
}
strtype newname = temp->Getname(*temp) + ROp.Getname(ROp) + ")";
temp->Nameit(newname);
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
VMatrix& Mlog(VMatrix &ROp)
{
// take logarithm of all elements
ROp.Garbage("Mlog");
VMatrix *temp = new VMatrix("log(", ROp.r, ROp.c);
for (int i = 1; i <= ROp.r; i++) {
for (int j = 1; j <= ROp.c; j++) {
double R = ROp.m(i, j);
temp->M(i, j) = (R > 0.0) ? log(R) :
Dispatch->Nrerror("Mlog: zero or negative argument to log");
}
}
strtype newname = temp->Getname(*temp) + ROp.Getname(ROp) + ")";
temp->Nameit(newname);
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
VMatrix& Msqrt(VMatrix &ROp)
{
// take sqrt of all elements
ROp.Garbage("Msqrt");
VMatrix *temp = new VMatrix("sqrt(", ROp.r, ROp.c);
for (int i = 1; i <= ROp.r; i++) {
for (int j = 1; j <= ROp.c; j++) {
double R = ROp.m(i, j);
temp->M(i, j) = (R >= 0.0) ? sqrt(R) :
Dispatch->Nrerror("Msqrt: zero or negative argument to sqrt");
}
}
strtype newname = temp->Getname(*temp) + ROp.Getname(ROp) + ")";
temp->Nameit(newname);
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
double Trace(VMatrix &ROp)
{
ROp.Garbage("Trace");
double trace = 0;
if (ROp.r != ROp.c)
ROp.Nrerror("Trace: matrix not square in trace");
for (int i = 1; i <= ROp.r; i++) trace += ROp.m(i, i);
return trace;
}
VMatrix& Ident(int n)
{
VMatrix *temp = new VMatrix("I(", n, n);
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) temp-> M(i, j) = (i == j) ? 1.0 : 0.0;
}
char buffer[MAXCHARS] = { '\0' };
gcvt(((double) n), NDECS, buffer);
strtype string = buffer;
string = temp->Getname(*temp) + string + ")";
temp->Nameit(string);
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
VMatrix& Helm(int n)
{
VMatrix *temp = new VMatrix("Helm(", n, n);
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
double d = 1.0 / sqrt((double) n);
double den = (j > 1) ? 1.0 / sqrt((double) j*(j - 1)) : d;
if (j > 1 && j < i) d = 0;
if (j > 1 && j == i) d = -den * ((double) (j - 1));
if (j > 1 && j > i) d = den;
temp->M(i, j) = d;
}
}
char buffer[MAXCHARS] = { '\0' };
gcvt(((double) n), NDECS, buffer);
strtype string = buffer;
string = temp->Getname(*temp) + string + ")";
temp->Nameit(string);
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
VMatrix& Fill(int r, int c, double d)
{
VMatrix *temp = new VMatrix("j", r, c);
for (int i = 1; i <= r; i++) {
for (int j = 1; j <= c; j++) temp-> M(i, j) = d;
}
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
VMatrix& Index(int lower, int upper, int step)
{
if (step == 0) Dispatch-> Nrerror("Index: step is zero");
int cnter = 0;
if (lower < upper) {
for (int i = lower; i <= upper; cnter++, i += step);
}
else {
for (int i = lower; i >= upper; cnter++, i += step);
}
if (cnter == 0)
Dispatch->Nrerror(
"Index: trying to allocate a matrix with zero rows");
VMatrix *temp = new VMatrix("Index", cnter, 1);
cnter = 1;
if (lower < upper) {
for (int i = lower; i <= upper; i += step)
temp->M(cnter++, 1) = (double) i;
}
else {
for (int i = lower; i >= upper; i += step)
temp->M(cnter++, 1) = (double) i;
}
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
VMatrix& Sweep(int k1, int k2, VMatrix& ROp) /* sweep out a row */
{
long double eps = ZERO, d;
int i, j, k, n, it;
ROp.Garbage("Sweep");
if (ROp.r != ROp.c) ROp.Nrerror("SWEEP: matrix a is not square");
n = ROp.r;
if (k2 < k1) { k = k1; k1 = k2; k2 = k; }
for (k = k1; k <= k2; k++) {
if (fabs(ROp.m(k, k)) < eps)
for (it = 1; it <= n; it++) {
ROp.M(it, k) = 0;
ROp.M(k, it) = 0;
}
else {
d = 1.0 / ROp.m(k, k);
ROp.M(k, k) = d;
for (i = 1; i <= n; i++) { if (i != k)
ROp.M(i, k) = ROp.m(i, k) *(double) - d; }
for (j = 1; j <= n; j++) { if (j != k)
ROp.M(k, j) = ROp.m(k, j) *(double) d; }
for (i = 1; i <= n; i++) {
if (i != k) {
for (j = 1; j <= n; j++) {
if (j != k)
ROp.M(i, j) = ROp.m(i, j) + (ROp.m(i, k)) *(ROp.m(k,j))/d;
}
}
}
}
}
strtype newname = "sweep(";
newname = newname + ROp.Getname(ROp) + ")";
ROp.Nameit(newname);
Dispatch->Push(ROp);
return Dispatch->ReturnMat();
} /*sweep*/
VMatrix& Inv(VMatrix &ROp) /*matrix inversion using sweep */
{
struct mat *b;
ROp.Garbage("Inv");
Dispatch->Inclevel();
if (ROp.c != ROp.r) ROp.Nrerror("INVERSE: matrix a not square");
strtype newname = "invs(";
newname = newname + ROp.Getname(ROp) + ")";
VMatrix *temp = new VMatrix("t", ROp.c, ROp.r);
*temp = ROp;
*temp = Sweep(1, temp->r, *temp);
temp->Nameit(newname);
Dispatch->Push(*temp);
delete temp;
return Dispatch->DecReturn();
} /* inverse */
VMatrix &Kron(VMatrix &a, VMatrix &b)
{
a.Garbage("Kron");
b.Garbage("Kron");
int cr = a.r*b.r;
int cc = a.c*b.c;
VMatrix *c = new VMatrix("(", cr, cc);
strtype mname = c->Getname(*c) + a.Getname(a) + "@" + b.Getname(b) + ")";
c->Nameit(mname);
for (int i = 1; i <= a.r; i++) {
for (int j = 1; j <= a.c; j++) {
for (int k = 1; k <= b.r; k++) {
for (int l = 1; l <= b.c; l++) {
c->M((i - 1) *b.r + k, (j - 1) *b.c + l) = a.m(i, j) *b.m(k,l);
}
}
}
}
Dispatch->Push(*c);
delete c;
return Dispatch->ReturnMat();
} /* kron */
void Pivot(VMatrix &Data, int RefRow,
double &Determ, enum boolean &DetEqualsZero)
{
DetEqualsZero = TTRUE;
int NewRow = RefRow;
while (DetEqualsZero && (NewRow < Data.r)) {
NewRow++;
if (fabs(Data.m(NewRow, RefRow)) > ZERO) {
for (int i = 1; i <= Data.r; i++) {
double temp = Data.m(NewRow, i);
Data.M(NewRow, i) = Data.m(RefRow, i);
Data.M(RefRow, i) = temp;
}
DetEqualsZero = FFALSE;
Determ = -Determ; /* row swap adjustment to det */
}
}
}
double Det(VMatrix &Datain)
{
Datain.Garbage("Det");
double Determ = 1;
if (Datain.r == Datain.c) {
Dispatch->Inclevel();
int Dimen = Datain.r;
VMatrix *Data = new VMatrix("data", Dimen, Dimen);
*Data = Datain;
long double Coef;
int Row, RefRow = 0;
enum boolean DetEqualsZero = FFALSE;
while (!(DetEqualsZero) && (RefRow < Dimen - 1)) {
RefRow++;
if (fabs(Data->m(RefRow, RefRow)) < ZERO)
Pivot(*Data, RefRow, Determ, DetEqualsZero);
if (!(DetEqualsZero))
for (Row = RefRow + 1; Row <= Dimen; Row++)
if (fabs(Data->m(Row, RefRow)) > ZERO) {
Coef = -((long double) Data->m(Row, RefRow)) /
((long double) Data->m(RefRow, RefRow));
for (int i = 1; i <= Dimen; i++)
Data->M(Row, i) = Data->m(Row, i) +
(double) (Coef*((long double) Data->m(RefRow,i)));
}
Determ *= Data->m(RefRow, RefRow);
}
Determ = (DetEqualsZero) ? 0 : Determ * Data-> m(Dimen, Dimen);
delete Data;
Dispatch->Declevel();
}
else Datain.Nrerror("Det: Matrix is not square\n");
return Determ;
}
//--------------- cholesky decomposition ----------------------
// ROp is symmetric, returns upper triangular S where ROp = S'S
//-------------------------------------------------------------
#define EPS 1.0e-7
VMatrix& Cholesky(VMatrix& ROp)
{
int nr = ROp.r;
ROp.Garbage("Cholesky");
for (int i = 1; i <= nr; i++) {
for (int j = 1; j < i; j++)
if (fabs(ROp.m(i, j) - ROp.m(j, i)) > EPS)
Dispatch->Nrerror("Cholesky: Input matrix is not symmetric");
}
VMatrix *temp = new VMatrix();
*temp = Fill(nr, nr, 0.0);
for (i = 1; i <= nr; i++) {
for (int j = i; j <= nr; j++) {
double sum = 0.0;
for (int k = 1; k < i; k++) {
sum += temp->m(k, i) * temp->m(k, j);
}
if (i == j) temp-> M(i, j) = (ROp.m(i, j) < sum) ?
Dispatch->Nrerror("Cholesky: negative pivot") :
sqrt(ROp.m(i, j) - sum);
else temp->M(i, j) = (fabs(temp->m(i, i)) < EPS) ?
Dispatch->Nrerror("Cholesky: zero or negative pivot") :
(ROp.m(i, j) - sum) / temp->m(i, i);
}
}
strtype newname = "half(";
newname = newname + ROp.Getname(ROp) + ")";
temp->Nameit(newname);
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
void QR(VMatrix& ROp, VMatrix& Q, VMatrix& R, boolean makeq)
{
ROp.Garbage("QR");
Q.Garbage("QR");
R.Garbage("QR");
Dispatch->Inclevel();
int r = ROp.r;
int c = ROp.c;
Q = ROp;
R = Fill(c, c, 0.0);
double s;
for (int j = 1; j <= c; j++) {
double sigma = 0.0;
for (int i = j; i <= r; i++) sigma += Q.m(i, j) *Q.m(i, j);
if (fabs(sigma) <= 1.0e-10) Dispatch->Nrerror("QR: singular X");
R.M(j, j) = s = (Q.m(j, j) < 0) ? sqrt(sigma) : - sqrt(sigma);
double beta = 1.0 /(s*Q.m(j, j) - sigma);
Q.M(j, j) = Q.m(j, j) - s;
for (int k = j + 1; k <= c; k++) {
double sum = 0.0;
for (int i = j; i <= r; i++) sum += Q.m(i, j) *Q.m(i, k);
sum *= beta;
for (i = j; i <= r; i++) Q.M(i, k) = Q.m(i, k) + Q.m(i, j) *sum;
R.M(j, k) = Q.m(j, k);
}
}
if (makeq) Q = ROp*Ginv(R);
Q.Nameit("Q");
R.Nameit("R");
Dispatch->Declevel();
}
#undef EPS
//------------------- Singular Valued Decomposition ----------
// from EISPACK SVD
//------------------------------------------------------------
long double safehypot(long double a1, long double b1)
// function to get hypotenuse numbers a1 and b1
// where a=log(a1) and b=log(b1)
{
long double del, sum;
long double a = 2.0*log(fabs(a1));
long double b = 2.0*log(fabs(b1));
del = (a > b) ? a : b;
// add a1^2 + b1^2, but in logarithms
sum = log(exp(a - del) + exp(b - del)) + del;
return (exp(0.5*sum));
}
void svdtemp(VMatrix &a, VMatrix &w, boolean matu, VMatrix &u,
boolean matv, VMatrix &v, int &ierr, VMatrix &rv1)
{
Dispatch->Inclevel();
a.Garbage("SVD");
w.Garbage("SVD");
u.Garbage("SVD");
v.Garbage("SVD");
rv1.Garbage("SVD");
int m = a.r;
int n = a.c;
int i, j, k, l, ii, i1, kk, k1, ll, l1, its, mn;
// VMatrix a("a",m,n),w("w",n,1),
// u("u",m,m),v("v",n,n),rv1("rv1",n,1);
long double c, f, g, h, s, x, y, z, eps, scale, machep, fss;
// boolean matu,matv;
// ********** Machep is a machine dependent parameter specifying
// the relative precision of floating point aritmetic.
//
// for pc doubles, range is 1.6^-308 to 1.6^308
// **************
machep = pow(1.6, -308.0);
ierr = 0;
u = a;
// ********householder reduction to bi diagonal form ********
g = 0.0;
scale = 0.0;
x = 0.0;
for (i = 1; i <= n; i++) {
l = i + 1;
rv1.M(i, 1) = scale*g;
g = 0.0;
s = 0.0;
scale = 0.0;
if (i > m) goto l210;
for (k = i; k <= m; k++) scale += fabs(u.m(k, i));
if (scale == 0.0) goto l210;
for (k = i; k <= m; k++) {
u.M(k, i) = u.m(k, i) / scale;
s += u.m(k, i) *u.m(k, i);
}
f = u.m(i, i);
g = -((f >= 0.0) ? fabs(sqrt(s)) : - fabs(sqrt(s)));
h = f*g - s;
u.M(i, i) = f - g;
if (i == n) goto l190;
for (j = l; j <= n; j++) {
s = 0.0;
for (k = i; k <= m; k++) s += u.m(k, i) *u.m(k, j);
f = s / h;
for (k = i; k <= m; k++) u.M(k, j) = u.m(k, j) + f*u.m(k, i);
}
l190 : for (k = i; k <= m; k++) u.M(k, i) = scale*u.m(k, i);
l210 : w.M(i, 1) = scale*g;
g = 0.0;
s = 0.0;
scale = 0.0;
if (i > m || i == n) goto l290;
for (k = l; k <= n; k++) scale += fabs(u.m(i, k));
if (scale == 0.0) goto l290;
for (k = l; k <= n; k++) {
u.M(i, k) = u.m(i, k) / scale;
s += u.m(i, k) *u.m(i, k);
}
f = u.m(i, l);
g = -((f >= 0.0) ? fabs(sqrt(s)) : - fabs(sqrt(s)));
h = f*g - s;
u.M(i, l) = f - g;
for (k = l; k <= n; k++) rv1.M(k, 1) = u.m(i, k) / h;
if (i == m) goto l270;
for (j = l; j <= m; j++) {
s = 0.0;
for (k = l; k <= n; k++) s += u.m(j, k) *u.m(i, k);
for (k = l; k <= n; k++) u.M(j, k) = u.m(j, k) + s*rv1.m(k, 1);
}
l270 : for (k = l; k <= n; k++) u.M(i, k) = scale*u.m(i, k);
l290 : x = (x > fabs(w.m(i, 1)) + fabs(rv1.m(i, 1))) ? x :
fabs(w.m(i, 1)) + fabs(rv1.m(i, 1));
}
// ********** accumulation of right-hand transformations ********
if (!matv) goto l410;
for (ii = 1; ii <= n; ii++) {
i = n + 1- ii;
if (i == n) goto l390;
if (g == 0.0) goto l360;
// double division avoids possible underflow
for (j = l; j <= n; j++) v.M(j, i) = (u.m(i, j) / u.m(i, l)) / g;
for (j = l; j <= n; j++) {
s = 0.0;
for (k = l; k <= n; k++) s += u.m(i, k) *v.m(k, j);
for (k = l; k <= n; k++) v.M(k, j) = v.m(k, j) + s*v.m(k, i);
}
l360 : for (j = l; j <= n; j++) {
v.M(i, j) = 0.0;
v.M(j, i) = 0.0;
}
l390 : v.M(i, i) = 1.0;
g = rv1.m(i, 1);
l = i;
}
// *************accumulation of left-hand transformations
l410 : if (!matu) goto l510;
// ************* for i= min(m,n) step -1 until 1 do ******
mn = n;
if (m < n) mn = m;
for (ii = 1; ii <= mn; ii++) {
i = mn + 1- ii;
l = i + 1;
g = w.m(i, 1);
if (i == n) goto l430;
for (j = l; j <= n; j++) u.M(i, j) = 0.0;
l430 : if (g == 0.0) goto l475;
if (i == mn) goto l460;
for (j = l; j <= n; j++) {
s = 0.0;
for (k = l; k <= m; k++) s += u.m(k, i) *u.m(k, j);
//double division avoids possible underflow
f =(s / u.m(i, i)) / g;
for (k = i; k <= m; k++) u.M(k, j) = u.m(k, j) + f*u.m(k, i);
}
l460 : for (j = i; j <= m; j++) u.M(j, i) = u.m(j, i) / g;
goto l490;
l475 : for (j = i; j <= m; j++) u.M(j, i) = 0.0;
l490 : u.M(i, i) = u.m(i, i) + 1.0;
}
// ********** diagonalization of the bidiagonal form ***********
l510 : eps = machep*x;
// ********** for k=n step -1 until 1 do ***********************
for (kk = 1; kk <= n; kk++) {
k1 = n - kk;
k = k1 + 1;
its = 0;
// ********** test for splitting .**********
// for l=k step -1 until 1 do
l520 : for (ll = 1; ll <= k; ll++) {
l1 = k - ll;
l = l1 + 1;
if (fabs(rv1.m(l, 1)) <= eps) goto l565;
// rv1(1) is always zero, so there is no exit
// through the bottom of the loop *********
if (fabs(w.m(l1, 1)) <= eps) goto l540;
}
// ************* cancellation of rv1(l) if l greater than 1******
l540 : c = 0.0;
s = 1.0;
for (i = l; i <= k; i++) {
f = s*rv1.m(i, 1);
rv1.M(i, 1) = c*rv1.m(i, 1);
if (fabs(f) <= eps) goto l565;
g = w.m(i, 1);
h = safehypot(f,g);
w.M(i, 1) = h;
c = g / h;
s = -f / h;
if (matu) { //go to 560
for (j = 1; j <= m; j++) {
y = u.m(j, l1);
z = u.m(j, i);
u.M(j, l1) = y*c + z*s;
u.M(j, i) = -y*s + z*c;
}
}
}
// ************** test for convergence ********************
l565 : z = w.m(k, 1);
if (l == k) goto l650;
// *************** if shift from bottom 2 by 2 minor *******
if (its == 50) goto l1000;
its = its + 1;
x = w.m(l, 1);
y = w.m(k1, 1);
g = rv1.m(k1, 1);
h = rv1.m(k, 1);
f =((y - z) *(y + z) + (g - h) *(g + h)) / (2.0*h*y);
g = safehypot(f,1.0);
fss =(f >= 0.0) ? fabs(g) : - fabs(g);
f =((x - z) *(x + z) + h*(y / (f + fss) - h)) / x;
// **************** next qr transformation ***************
c = 1.0;
s = 1.0;
//
for (i1 = l; i1 <= k1; i1++) {
i = i1 + 1;
g = rv1.m(i, 1);
y = w.m(i, 1);
h = s*g;
g = c*g;
z = safehypot(f,h);
rv1.M(i1, 1) = z;
c = f / z;
s = h / z;
f = x*c + g*s;
g = -x*s + g*c;
h = y*s;
y = y*c;
if (matv) { // goto l575;
for (j = 1; j <= n; j++) {
x = v.m(j, i1);
z = v.m(j, i);
v.M(j, i1) = x*c + z*s;
v.M(j, i) = -x*s + z*c;
}
}
z = safehypot(f,h);
w.M(i1, 1) = z;
// ************ rotation can be arbitrary if z is zero *******
if (z != 0.0) { // goto l580;
c = f / z;
s = h / z;
}
// l580:
f = c*g + s*y;
x = -s*g + c*y;
if (matu) { //goto l600;
for (j = 1; j <= m; j++) {
y = u.m(j, i1);
z = u.m(j, i);
u.M(j, i1) = y*c + z*s;
u.M(j, i) = -y*s + z*c;
}
}
l600 : x = x; } //continue
rv1.M(l, 1) = 0.0;
rv1.M(k, 1) = f;
w.M(k, 1) = x;
goto l520;
// ************** convergence **************
l650 : if (z >= 0.0) goto l700;
// ************* w(k) is made non-negative *********
w.M(k, 1) = -z;
if (!matv) goto l700;
for (j = 1; j <= n; j++) v.M(j, k) = -v.m(j, k);
l700 : x = x; }
ierr = 0;
Dispatch->Declevel();
return;
// ***************** set error -- no convergence to a
// singular value after 30 iterations
l1000 : ierr = k;
Dispatch->Declevel();
return;
}
int Svd(VMatrix &A, VMatrix &U, VMatrix &S, VMatrix &V,
boolean makeu, boolean makev)
{
Dispatch->Inclevel();
A.Garbage("SVD");
VMatrix aa, uu, ss, vv, rv1;
int ierr = 0;
if (A.r < A.c) aa = Tran(A);
else aa = A;
uu = Fill(aa.r, aa.r, 0.0);
vv = Fill(aa.c, aa.c, 0.0);
ss = Fill(aa.c, 1, 0.0);
rv1 = Fill(aa.c, 1, 0.0);
svdtemp(aa, ss, makeu, uu, makev, vv, ierr, rv1);
if (A.r < A.c) {
if (makeu) U = vv;
if (makev) V = uu;
} else {
if (makeu) U = uu;
if (makev) V = vv;
}
S = ss;
Dispatch->Declevel();
return ierr;
}
//---------------- end svd ------------------------------------
VMatrix& Ginv(VMatrix& ROp)
{
ROp.Garbage("Ginv");
Dispatch->Inclevel();
VMatrix *u = new VMatrix;
VMatrix *s = new VMatrix;
VMatrix *v = new VMatrix;
VMatrix *g = new VMatrix;
Svd(ROp, *u, *s, *v);
for (int i = 1; i <= s-> r; i++) {
double t = s->m(i, 1);
s->M(i, 1) = (fabs(t) <= 0.0) ? 0.0 : 1.0 / t;
}
*g = *v *Diag(*s) *Tran(*u);
strtype newname = "ginv(";
newname = newname + ROp.Getname(ROp) + ")";
g->Nameit(newname);
Dispatch->Push(*g);
delete u;
delete s;
delete v;
delete g;
return Dispatch->DecReturn();
}
//-------------------------- fft ------------------------------
// de Boor (1980) SIAM J SCI. STAT. COMPUT. V1 no 1, pp 173-178
// and NewMat by R. G. Davies a C++ matrix Package
//-------------------------------------------------------------
static void cossin(int n, int d, double& c, double& s)
// calculate cos(twopi*n/d) and sin(twopi*n/d)
// minimise roundoff error
{
long n4 = n * 4;
int sector =(int) floor((double) n4 / (double) d + 0.5);
n4 -= sector * d;
if (sector < 0) sector = 3 - (3 - sector) % 4; else sector %= 4;
double ratio = 1.5707963267948966192 * (double) n4 / (double) d;
switch (sector) {
case 0 : c = cos(ratio);
s = sin(ratio);
break;
case 1 : c = -sin(ratio);
s = cos(ratio);
break;
case 2 : c = -cos(ratio);
s = -sin(ratio);
break;
case 3 : c = sin(ratio);
s = -cos(ratio);
break;
}
}
static void fftstep(VMatrix &AB, VMatrix &XY,
int after, int now, int before)
{
int gamma = after * before;
int delta = now * after;
double r_value, i_value, ridx_value, iidx_value, temp;
int j, x1, ia, a, a1, x2, ib, a2, idx, in;
double r_arg = 1.0, i_arg = 0.0;
int x = 1;
int m = AB.r - gamma;
for (j = 0; j < now; j++) {
a = 1;
x1 = x;
x += after;
for ( ia = 0; ia < after; ia++) {
// generate sins & cosines explicitly rather than iteratively
// for more accuracy; but slower
cossin(- (j*after + ia), delta, r_arg, i_arg);
a1 = a++;
x2 = x1++;
if (now == 2) {
ib = before;
while (ib--) {
a2 = m + a1;
a1 += after;
idx = a2 - gamma;
r_value = AB.m(a2, 1);
i_value = AB.m(a2, 2);
ridx_value = AB.m(idx,1);
iidx_value = AB.m(idx,2);
XY.M(x2, 1) = r_value*r_arg - i_value*i_arg + ridx_value;
XY.M(x2, 2) = r_value*i_arg + i_value*r_arg + iidx_value;
x2 += delta;
}
}
else {
ib = before;
while (ib--) {
a2 = m + a1;
a1 += after;
r_value = AB.m(a2, 1);
i_value = AB.m(a2, 2);
in = now - 1;
while (in--) {
a2 -= gamma;
ridx_value = AB.m(a2, 1);
iidx_value = AB.m(a2, 2);
temp = r_value;
r_value = r_value*r_arg - i_value*i_arg + ridx_value;
i_value = temp*i_arg + i_value*r_arg + iidx_value;
}
XY.M(x2, 1) = r_value;
XY.M(x2, 2) = i_value;
x2 += delta;
}
}
}
}
}
VMatrix& Fft(VMatrix &ROp, boolean INZEE)
// if INZEE == TTRUE then calculate fft
// if INZEE == FFALSE then calculate inverse fft
{
ROp.Garbage("FFT");
if (ROp.c > 2) Dispatch->Nrerror("FFT: Input has more than two columns");
int n = ROp.r; // length of arrays
const int nextmx = 26;
int prime[26] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37,
41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83,
89, 97, 101 };
int after = 1;
int before = n;
int next = 0;
boolean inzee = TTRUE;
Dispatch->Inclevel();
VMatrix *AB = new VMatrix("AB", n, 2);
VMatrix *XY = new VMatrix("XY", n, 2);
if (ROp.c == 1) {
for (int i = 1; i <= n; i++) {
AB->M(i, 1) = ROp.m(i, 1);
AB->M(i, 2) = 0.0;
}
} else *AB = ROp;
*XY = Fill(n, 2, 0.0);
if (INZEE == FFALSE) {
// take complex conjugate for ifft
for (int i = 1; i <= n; i++) AB-> M(i, 2) = -AB->m(i, 2);
}
do {
int now, b1;
for (;;) {
if (next < nextmx) now = prime[next];
b1 = before / now;
if (b1 * now == before) break;
next++;
now += 2;
}
before = b1;
if (inzee) fftstep(*AB, *XY, after, now, before);
else fftstep(*XY, *AB, after, now, before);
inzee =(inzee == TTRUE) ? FFALSE : TTRUE;
after *= now;
}
while (before != 1);
if (inzee == TTRUE) *XY = *AB;
// divide by n for ifft
if (INZEE == FFALSE) *XY = (*XY) / ((double) n);
strtype newname = (INZEE == TTRUE) ? "fft( " : "ifft(";
newname = newname + ROp.Getname(ROp) + ")";
XY->Nameit(newname);
Dispatch->Push(*XY);
delete XY;
delete AB;
return Dispatch->DecReturn();
}
//////////////////// reshaping functions
VMatrix& Vec(VMatrix& ROp)
{ // send columns of ROp to a vector
ROp.Garbage("Vec");
long lrc = ((long) ROp.r) *((long) ROp.c);
if (lrc >= 32768 || lrc < 1)
Dispatch->Nrerror("Vec: vec produces invalid indices");
int r = ROp.r*ROp.c;
int c = 1;
VMatrix *temp = new VMatrix("t", r, c);
int cnter = 1;
for (int j = 1; j <= ROp.c; j++)
for (int i = 1; i <= ROp.r; i++) temp->M(cnter++, 1) = ROp.m(i, j);
strtype newname = "vec(";
newname = newname + ROp.Getname(ROp) + ")";
temp->Nameit(newname);
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
VMatrix& Vecdiag(VMatrix& ROp)
{ // Extract the diagonal from a square matrix and put in a vector
ROp.Garbage("Vecdiag");
if (ROp.r != ROp.c)
Dispatch->Nrerror("Vecdiag: ROp is not square");
VMatrix *temp = new VMatrix("vdiag(", ROp.r, 1);
for (int i = 1; i <= ROp.r; i++) temp->M(i, 1) = ROp.m(i, i);
strtype newname = "vdiag(";
newname = newname + ROp.Getname(ROp) + ")";
temp->Nameit(newname);
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
VMatrix& Diag(VMatrix& ROp)
{ // make a diagonal matrix from a vector or the principal diagonal of
// a square matrix, off diagonal elements are zero
ROp.Garbage("Diag");
if (ROp.r != ROp.c && ROp.c != 1)
Dispatch->Nrerror(
"Diag: ROp is not square or is not a column vector" );
Dispatch->Inclevel();
VMatrix *temp = new VMatrix;
*temp = Fill(ROp.r, ROp.r, 0.0);
int rr = ROp.r;
int cc = ROp.c;
for (int i = 1; i <= rr; i++)
temp->M(i, i) = (rr == cc) ? ROp.m(i, i) : ROp.m(i, 1);
strtype newname = "diag(";
newname = newname + ROp.Getname(ROp) + ")";
temp->Nameit(newname);
Dispatch->Push(*temp);
delete temp;
return Dispatch->DecReturn();
}
VMatrix& Shape(VMatrix& ROp, int nrows)
{ // reshape a matrix into n rows, nrows must divide r*c without a
// remainder
ROp.Garbage("Shape");
long nr = (long) nrows;
long lr = (long) ROp.r;
long lc = (long) ROp.c;
if (lr*lc % nr)
Dispatch->Nrerror("Shape: nrows divides r*c with a remainder");
Dispatch->Inclevel();
VMatrix *temp = new VMatrix;
*temp = ROp;
long lrc = (lr*lc) / nr;
if (nr >= 32768 || nr < 1 || lrc >= 32768 || lrc < 1)
Dispatch->Nrerror("Shape: reshape produces invalid indices");
temp->r = nrows;
temp->c = (int) lrc;
strtype newname = "shape(";
newname = newname + ROp.Getname(ROp) + ")";
temp->Nameit(newname);
Dispatch->Push(*temp);
delete temp;
return Dispatch->DecReturn();
}
///////////////////////////// summation functions // maxs and mins
//// typedef enum margins { ALL, ROWS, COLUMNS } margins;
VMatrix& Sum(VMatrix& ROp, margins marg)
{ // sum(ROp, ALL) returns 1x1
// sum(ROp, ROWS) returns 1xc
// sum(ROp, COLUMNS) returns cx1
ROp.Garbage("Sum");
int r, c;
double sum;
switch (marg) {
case ALL : r = 1; c = 1; break;
case ROWS : r = 1; c = ROp.c; break;
case COLUMNS : r = ROp.r; c = 1; break;
default : Dispatch->Nrerror("Sum: invalid margin specification");
}
VMatrix *temp = new VMatrix("t", r, c);
int i, j;
switch (marg) {
case ALL : sum = 0.0;
for (i = 1; i <= ROp.r; i++)
for (j = 1; j <= ROp.c; j++) sum += ROp.m(i, j);
temp->M(1, 1) = sum;
break;
case ROWS : for (j = 1; j <= c; j++) {
sum = 0.0;
for (i = 1; i <= ROp.r; i++) sum += ROp.m(i, j);
temp->M(1, j) = sum;
}
break;
case COLUMNS : for (i = 1; i <= r; i++) {
sum = 0.0;
for (j = 1; j <= ROp.c; j++) sum += ROp.m(i, j);
temp->M(i, 1) = sum;
}
break;
}
strtype newname = "sum(";
newname = newname + ROp.Getname(ROp) + ")";
temp->Nameit(newname);
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
VMatrix& Sumsq(VMatrix& ROp, margins marg)
{ // sumsq(ROp, ALL) returns 1x1
// sumsq(ROp, ROWS) returns 1xc
// sumsq(ROp, COLUMNS) returns cx1
ROp.Garbage("Sum");
int r, c;
double sum;
switch (marg) {
case ALL : r = 1; c = 1; break;
case ROWS : r = 1; c = ROp.c; break;
case COLUMNS : r = ROp.r; c = 1; break;
default : Dispatch->Nrerror("Sumsq: invalid margin specification");
}
VMatrix *temp = new VMatrix("t", r, c);
int i, j;
switch (marg) {
case ALL : sum = 0.0;
for (i = 1; i <= ROp.r; i++)
for (j = 1; j <= ROp.c; j++) sum += ROp.m(i, j) *ROp.m(i, j);
temp->M(1, 1) = sum;
break;
case ROWS : for (j = 1; j <= c; j++) {
sum = 0.0;
for (i = 1; i <= ROp.r; i++) sum += ROp.m(i, j) *ROp.m(i, j);
temp->M(1, j) = sum;
}
break;
case COLUMNS : for (i = 1; i <= r; i++) {
sum = 0.0;
for (j = 1; j <= ROp.c; j++) sum += ROp.m(i, j) *ROp.m(i, j);
temp->M(i, 1) = sum;
}
break;
}
strtype newname = "sumsq(";
newname = newname + ROp.Getname(ROp) + ")";
temp->Nameit(newname);
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
VMatrix& Cusum(VMatrix& ROp)
{ // cumulative sum along rows
ROp.Garbage("Cusum");
VMatrix *temp = new VMatrix("t", ROp.r, ROp.c);
double sum = 0.0;
for (int i = 1; i <= ROp.r; i++) {
for (int j = 1; j <= ROp.c; j++) {
sum += ROp.m(i, j);
temp->M(i, j) = sum;
}
}
strtype newname = "cusum(";
newname = newname + ROp.Getname(ROp) + ")";
temp->Nameit(newname);
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
VMatrix& Mmax(VMatrix& ROp, margins marg)
{ // Mmax(ROp, ALL) returns 3x1
// Mmax(ROp, ROWS) returns 3xc
// Mmax(ROp, COLUMNS) returns cx3
ROp.Garbage("Sum");
int r, c;
switch (marg) {
case ALL : r = 3; c = 1; break;
case ROWS : r = 3; c = ROp.c; break;
case COLUMNS : r = ROp.r; c = 3; break;
default : Dispatch->Nrerror("Mmax: invalid margin specification");
}
VMatrix *temp = new VMatrix("t", r, c);
int maxr, maxc, i, j;
double mmax;
switch (marg) {
case ALL : mmax = ROp.m(1, 1);
maxr = 1; maxc = 1;
for (i = 1; i <= ROp.r; i++)
for (j = 1; j <= ROp.c; j++)
mmax = (mmax < ROp.m(i, j)) ?
maxr = i, maxc = j, ROp.m(i, j) : mmax;
temp->M(1, 1) = (double) maxr;
temp->M(2, 1) = (double) maxc;
temp->M(3, 1) = mmax;
break;
case ROWS : for (j = 1; j <= c; j++) {
mmax = ROp.m(1, j);
maxr = 1; maxc = j;
for (i = 1; i <= ROp.r; i++)
mmax = (mmax < ROp.m(i, j)) ? maxr = i, ROp.m(i,j): mmax;
temp->M(1, j) = (double) maxr;
temp->M(2, j) = (double) maxc;
temp->M(3, j) = mmax;
}
break;
case COLUMNS : for (i = 1; i <= r; i++) {
mmax = ROp.m(i, 1);
maxr = i; maxc = 1;
for (j = 1; j <= ROp.c; j++)
mmax = (mmax < ROp.m(i, j)) ? maxc = j, ROp.m(i,j): mmax;
temp->M(i, 1) = (double) maxr;
temp->M(i, 2) = (double) maxc;
temp->M(i, 3) = mmax;
}
break;
}
strtype newname = "max(";
newname = newname + ROp.Getname(ROp) + ")";
temp->Nameit(newname);
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
VMatrix& Mmin(VMatrix& ROp, margins marg)
{ // Mmin(ROp, ALL) returns 3x1
// Mmin(ROp, ROWS) returns 3xc
// Mmin(ROp, COLUMNS) returns cx3
ROp.Garbage("Sum");
int r, c;
switch (marg) {
case ALL : r = 3; c = 1; break;
case ROWS : r = 3; c = ROp.c; break;
case COLUMNS : r = ROp.r; c = 3; break;
default : Dispatch->Nrerror("Mmin: invalid margin specification");
}
VMatrix *temp = new VMatrix("t", r, c);
int maxr, maxc, i, j;
double mmin;
switch (marg) {
case ALL : mmin = ROp.m(1, 1);
maxr = 1; maxc = 1;
for (i = 1; i <= ROp.r; i++)
for (j = 1; j <= ROp.c; j++)
mmin = (mmin > ROp.m(i, j)) ?
maxr = i, maxc = j, ROp.m(i, j) : mmin;
temp->M(1, 1) = (double) maxr;
temp->M(2, 1) = (double) maxc;
temp->M(3, 1) = mmin;
break;
case ROWS : for (j = 1; j <= c; j++) {
mmin = ROp.m(1, j);
maxr = 1; maxc = j;
for (i = 1; i <= ROp.r; i++)
mmin = (mmin > ROp.m(i, j)) ? maxr = i, ROp.m(i,j): mmin;
temp->M(1, j) = (double) maxr;
temp->M(2, j) = (double) maxc;
temp->M(3, j) = mmin;
}
break;
case COLUMNS : for (i = 1; i <= r; i++) {
mmin = ROp.m(i, 1);
maxr = i; maxc = 1;
for (j = 1; j <= ROp.c; j++)
mmin = (mmin > ROp.m(i, j)) ? maxc = j, ROp.m(i,j): mmin;
temp->M(i, 1) = (double) maxr;
temp->M(i, 2) = (double) maxc;
temp->M(i, 3) = mmin;
}
break;
}
strtype newname = "min(";
newname = newname + ROp.Getname(ROp) + ")";
temp->Nameit(newname);
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
//////////////////////////////// elemenatary row and column operations
////////////////// note these functions operate directly on ROp
void Crow(VMatrix& ROp, int row, double c)
{ // multiply a row by a constant
ROp.Garbage("Crow");
if (row < 1 || row > ROp.r)
Dispatch->Nrerror("Crow: row index out of range");
for (int i = 1; i <= ROp.c; i++) ROp.M(row, i) = c*ROp.m(row, i);
return;
}
void Srow(VMatrix &ROp, int row1, int row2)
{ // swap rows
ROp.Garbage("Srow");
if (row1 < 1 || row1 > ROp.r || row2 < 1 || row2 > ROp.r)
Dispatch->Nrerror("Srow: row index out of range");
double tmp = 0;
for (int i = 1; i <= ROp.c; i++) {
tmp = ROp.m(row1, i);
ROp.M(row1, i) = ROp.m(row2, i);
ROp.M(row2, i) = tmp;
}
return;
}
void Lrow(VMatrix &ROp, int row1, int row2, double c)
{ // add c*r1 to r2
ROp.Garbage("Srow");
if (row1 < 1 || row1 > ROp.r || row2 < 1 || row2 > ROp.r)
Dispatch->Nrerror("Lrow: row index out of range");
for (int i = 1; i <= ROp.c; i++)
ROp.M(row2, i) = ROp.m(row2, i) + c*ROp.m(row1, i);
return;
}
void Ccol(VMatrix& ROp, int col, double c)
{ // multiply a col by a constant
ROp.Garbage("Ccol");
if (col < 1 || col > ROp.c)
Dispatch->Nrerror("Ccol: col index out of range");
for (int i = 1; i <= ROp.r; i++) ROp.M(i, col) = c*ROp.m(i, col);
return;
}
void Scol(VMatrix &ROp, int col1, int col2)
{ // swap cols
ROp.Garbage("Scol");
if (col1 < 1 || col1 > ROp.c || col2 < 1 || col2 > ROp.c)
Dispatch->Nrerror("Scol: col index out of range");
double tmp = 0;
for (int i = 1; i <= ROp.r; i++) {
tmp = ROp.m(i, col1);
ROp.M(i, col1) = ROp.m(i, col2);
ROp.M(i, col2) = tmp;
}
return;
}
void Lcol(VMatrix &ROp, int col1, int col2, double c)
{ // add c*c1 to c2
ROp.Garbage("Scol");
if (col1 < 1 || col1 > ROp.c || col2 < 1 || col2 > ROp.c)
Dispatch->Nrerror("Lcol: col index out of range");
for (int i = 1; i <= ROp.r; i++)
ROp.M(i, col2) = ROp.m(i, col2) + c*ROp.m(i, col1);
return;
}